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Abstract — Model-based prognostics approaches employ do- 
main knowledge about a system, its components, and how 
they fail through the use of physics-based models. Compo- 
nent wear is driven by several different degradation phenom- 
ena, each resulting in their own damage progression path, 
overlapping to contribute to the overall degradation of the 
component. We develop a model-based prognostics method- 
ology using particle filters, in which the problem of charac- 
terizing multiple damage progression paths is cast as a joint 
state-parameter estimation problem. The estimate is repre- 
sented as a probability distribution, allowing the prediction 
of end of life and remaining useful life within a probabilistic 
framework that supports uncertainty management. We also 
develop a novel variance control mechanism that maintains an 
uncertainty bound around the hidden parameters to limit the 
amount of estimation uncertainty and, consequently, reduce 
prediction uncertainty. We construct a detailed physics-based 
model of a centrifugal pump, to which we apply our model- 
based prognostics algorithms. We illustrate the operation of 
the prognostic solution with a number of simulation-based 
experiments and demonstrate the performance of the chosen 
approach when multiple damage mechanisms are active. 
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1. Introduction 

Model-based prognostics approaches employ domain knowl- 
edge about a system, its components, and how they fail 
through the use of physics-based models that capture the 
underlying physical phenomena [1-3]. Component wear is 
driven by several different degradation phenomena. Each of 
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these degradation phenomena results in its own damage pro- 
gression path, which all combine to contribute to the over- 
all degradation of the component. Due to manufacturing 
variances and differences in usage and environmental con- 
ditions, the damage progression rates for the different dam- 
age mechanisms vary among components of the same type. 
This poses considerable challenges to data-driven (model- 
free) approaches, which use run-to-failure data to train ma- 
chine learning algorithms to make end of life and remaining 
useful life predictions [4], because often the training data to 
cover a sufficient portion of such cases is lacking. In the ab- 
sence of such data, model-based approaches are better-suited, 
since they use underlying physical models to help estimate 
the amount of damage and the rates of damage progression. 

Extending previous work in [1], we develop a model-based 
prognostics methodology using particle filters, in which the 
problem of characterizing multiple damage progression paths 
is cast as a joint state-parameter estimation problem. The 
estimate is represented as a probability distribution, allow- 
ing the prediction of end of life and remaining useful life 
within a probabilistic framework that supports uncertainty 
management. In particle filter-based parameter estimation, 
an artificial random walk evolution is assigned to the parame- 
ters, which is necessary for convergence of the estimates and 
proper tracking afterwards. But, the optimal variance of the 
random walk depends on the unknown parameter value. To 
reduce the amount of this artificial uncertainty, we introduce 
a novel variance control mechanism that maintains an uncer- 
tainty bound around an unknown parameter being estimated. 

We demonstrate our prognostics methodology on a centrifu- 
gal pump. Centrifugal pumps are used in a wide range of 
applications, from water supply to spacecraft fueling systems. 
Because pumps typically see high usage, they can particularly 
benefit from prognostics and health management solutions to 
ensure satisfactory system performance, extended component 
lifetime, and limited downtime. Model-based diagnosis has 
been investigated previously with centrifugal pumps [5-7]. 
However, most prognostics approaches for pumps have been 
data-driven, usually based on pump vibration signals. A 
principal component analysis method is applied for condi- 
tion monitoring of a pump using vibration signals in [8]. A 
model-based approach is presented in [9], however it consid- 
ers only a single degradation mode. We illustrate here our 



model-based prognostic approach for centrifugal pumps us- 
ing a number of simulation-based experiments when multiple 
damage mechanisms are active. We evaluate algorithm per- 
formance using established prognostics metrics [10]. 

The paper is organized as follows. Section 2 formally defines 
the prognostics problem and describes the prognostics archi- 
tecture. Section 3 describes the modeling methodology and 
develops the centrifugal pump model for prognostics. Sec- 
tion 4 describes the particle filter-based damage estimation 
method and develops the variance control scheme. Section 5 
discusses the prediction methodology. Section 6 provides 
results from a number of simulation-based experiments and 
evaluates the approach. Section 7 concludes the paper. 

2. Prognostics Approach 

The problem of prognostics is to predict the EOT and/or the 
RUL of a component. In this section, we first formally de- 
fine the problem of prognostics. We then describe a general 
model-based architecture for prognostics. 

Problem Formulation 

In general, a system model may be defined as 

x(t) = f(t,x(/),0(/),u(/),v(/)) 

y(t) = h(f, x(t), 0(f), u(/), n(/)), 

where x(/) G is the state vector, 6{t) G M"® is the 
parameter vector, u(f) G M"“ is the input vector, v(t) G K"’' 
is the process noise vector, f is the state equation, y{t) G M"* 
is the output vector, n(f) G M"" is the measurement noise 
vector, and h is the output equation. This form represents a 
general nonlinear model with no restrictions on the functional 
forms of f or h. Further, the noise terms may be coupled in a 
nonlinear way with the states and parameters. The parameters 
9{t) evolve in an unknown way, but are typically considered 
to be constant in practice. 

The goal is to predict EOL (and/or RUL) at a given time 
point tp using the discrete sequence of observations up to 
time tp, denoted as y^-tp- EOL is defined as the time point 
at which the component no longer meets a functional re- 
quirement (e.g., a pump is overheated). This point is often 
linked to a damage threshold, beyond which the component 
fails to function properly. In general, we may express this 
threshold as a function of the system state and parameters, 
TpoLi^it), 9{t)), which determines whether EOL has been 
reached, where 

Tbol(x(/),0(/)) = | J’ 


if EOL is reached 
otherwise. 


The EOL threshold is linked to a boundary in the 
multi-dimensional damage space. Inside the bound- 
7bol(x(/), 0(f)) = 0, and outside the boundary, 

TEOL{^{t),d{t)) = 1. Fig. 1 illustrates this concept with 
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a two-dimensional example, with damage dimensions di and 
d 2 - The dimensions are normalized such that d\ = \ cor- 
responds to the maximum allowable damage for d\ when 
d 2 = 0, and ^2 = 1 corresponds to the maximum allowable 
damage for c ?2 when d\ = 0. If the different damage mech- 
anisms are considered independently, then the space where 
7bol(x(/), 9{t)) = 0 would be defined by the space within 
the dashed lines in the figure. In higher dimensions, this 
space forms a hypercube. However, in general, the different 
damage mechanisms cannot be considered independently in 
defining EOL, because increased damage along one dimen- 
sion may either allow a greater amount of damage or restrict 
the allowable amount of damage along another damage di- 
mension. For example, in a normally-closed valve, where 
EOL is defined by opening and closing times, friction damage 
will cause the valve to open more slowly, but a weakening of 
the return spring will allow the valve to open more quickly. 
So, the actual EOL threshold may take on a more complex 
form, as shown by the shaded area in Fig. 1 . In the regions 
of the space where Tj;ol(x(/), 9{t)) = 0 that extend beyond 
the hypercube, more damage is allowed, and in the regions 
that fall within the hypercube, damage is restricted further. 

Using Teol, we can formally define EOL with 

EOL{tp) = argminTEOL(x(/),0(/)) = 1, 

t>tp 

i.e., EOL is the earliest time point at which the damage 
threshold is met. RUL may then be defined with 

RUL(tp) = EOL{tp)-tp. 

Note that we are interested in the EOL formed by the com- 
bined effects of all damage progressions paths, so they must 
be considered simultaneously, rather than independently. 

In practice, many sources of uncertainty exist that affect the 
prediction. Noise is inherent in the process and the mea- 
surements, represented by the noise terms v(/) and n(f), re- 
spectively. Further, the future inputs of the system, which 
affect the evolution of the state, and therefore the progres- 
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Figure 2. Prognostics architecture. 


sion of damage, are not always known. Certain input pro- 
files may also excite some damage mechanisms more than 
others. Thus, it is much more useful to compute a proba- 
bility distribution of the EOL or RUL, rather than a single 
prediction point. The goal, then, is to compute, at time tp, 

p{EOL{tp)\yo,tj.) or p{RUL{tp)\yo:tp)- 

Prognostics Architecture 

In our model-based approach, we develop detailed physics- 
based models of components and systems that include de- 
scriptions of how fault parameters evolve in time. These 
models depend on unknown and possibly time-varying wear 
parameters, 9{t). Therefore, our solution to the prognostics 
problem takes the perspective of joint state-parameter estima- 
tion. In discrete time k, we estimate and 9k, and use these 
estimates to predict EOL and RUL at desired time points. 

We employ the prognostics architecture in Eig. 2. The sys- 
tem is provided with inputs Uj. and provides measured out- 
puts yfc. Prognostics may begin at t = 0, with the dam- 
age estimation module determining estimates of the states 
and unknown parameters, represented as a probability dis- 
tribution p{xk,9k\yo:k)- In parallel, a fault detection, iso- 
lation, and identification (EDII) module may be used to de- 
termine which damage mechanisms are active, represented as 
a fault set F. The damage estimation module may then use 
this result to limit the space of parameters that must be esti- 
mated. Alternatively, prognostics may begin only when diag- 
nostics has completed. The prediction module uses the joint 
state-parameter distribution, along with hypothesized future 
inputs, to compute EOL and RUL as probability distributions 
p{EOLkp\yo-.kp) and p{RULkp\yo-.kp) at given prediction 
times kp. In this paper, we focus on the damage estimation 
and prediction modules, and assume that the EDII module 
does not inform the prognostics, i.e., all possible damage pro- 
gression paths must be tracked starting from t = 0. 

3. Pump Modeling 

We apply our prognostics approach to a centrifugal pump, 
and develop a physics-based model of its nominal and faulty 
behavior. Centrifugal pumps are used in a variety of do- 
mains for fluid delivery. A schematic of a typical centrifugal 
pump is shown in Eig. 3. Eluid enters the inlet, and the rota- 
tion of the impeller forces fluid through the outlet. The im- 
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peller is driven by an electric motor, typically a three-phase 
alternating-current induction motor. The radial and thrust 
bearings help to minimize friction along the pump shaft. The 
bearing housing contains oil which lubricates the bearings. A 
seal prevents fluid flow into the bearing housing. Wear rings 
prevent internal pump leakage from the outlet to the inlet side 
of the impeller, but a small clearance is typically allowed to 
minimize friction (a small internal leakage is normal). 

The state of the pump is given by 

x(f)=[u;(f) Tt(t) Tr{t) To(t)f , 

where uj{t) is the rotational velocity of the pump, Tt{t) is the 
thrust bearing temperature, (f) is the radial bearing temper- 
ature, and To(f) is the oil temperature. 

The rotational velocity of the pump is described using a 
torque balance, 

W = j (Te(f) - ruj{t) - TL{t)) , 

where J is the lumped motor/pump inertia, Tg is the electro- 
magnetic torque provided by the motor, r is the lumped fric- 
tion parameter, and Tp is the load torque. In an induction mo- 
tor, a voltage is applied to the stationary part, the stator, which 
creates a current through the stator coils. With a polyphase 
supply, this creates a rotating magnetic field which induces 
a current in the rotating part, the rotor, causing it to turn. A 
torque is produced on the rotor only when there is a differ- 
ence between the synchronous speed of the supply voltage, 
LOs and the mechanical rotation, uj. This difference, called 
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Figure 4. Induction motor equivalent circuit. 


slip, is defined as 

UJs — U) 

s = . 

UJs 

The expression for the torque Te is derived from an equiva- 
lent circuit representation for the three-phase induction mo- 
tor, shown in Fig. 4, based on rotor and stator resistances and 
inductances, and the slip s [11]: 

^rms 

su>s {Ri + R 2 / + {uJsLi uigL2)‘^ 

where i?i is the stator resistance, Li is the stator inductance, 
i ?2 is the rotor resistance, L 2 is the rotor inductance, n is 
the number of phases (typically 3), and p is the number of 
magnetic pole pairs. For a 3600 rpm motor, p = 1. The de- 
pendence of torque on slip creates a feedback loop that causes 
the rotor to follow the rotation of the magnetic field. The ro- 
tor speed may be controlled by changing the input frequency 
UJs, e.g., through the use of a variable-frequency drive. 

The load torque is a polynomial function of the flow rate 
through the pump and the impeller rotational velocity [5,6]: 

tl = a^uo^ + aiUjQ - a2Q^, 

where Q is the flow, and qq, ai, and 02 are coefficients de- 
rived from the pump geometry [6] . 

The rotation of the impeller creates a pressure difference from 
the inlet to the outlet of the pump, which drives the pump 
flow, Q. The pump pressure is computed as 

Pp = Auj^ -f b\ujQ — & 2 Q ^7 

where A is the impeller area, and b\ and &2 are coefficients 
derived from the pump geometry. Flow through the impeller, 
Qi, is computed using the pressure differences: 

Qi = c^\Ps +Pp -pd\sign{ps + Pp - Pd), 

where c is a flow coefficient, ps is the suction pressure, and 
Pd is the discharge pressure. The small (normal) leakage flow 
from the discharge end to the suction end due to the clearance 
between the wear rings and the impeller is described by 

Qi = cis/\pd- Ps\sign{pd -Ps), 


where ci is a flow coefficient. The discharge flow, Q, is then 

Q — Qi Qi ■ 

Pump temperatures are often monitored as indicators of pump 
condition. The oil heats up due to the radial and thrust bear- 
ings and cools to the environment: 

To = ^{Ho,l{Tt - To) + HoMTr - To) - Ho^To ~ Ta)), 
'Jo 

where Jo is the thermal inertia of the oil, and the Ho^i terms 
are heat transfer coefficients. The thrust bearings heat up due 
to the friction between the pump shaft and the bearings, and 
cool to the oil and the environment: 

ft = - Ht,i{Tt - To) - Ht,2{Tt - T„)), 

•Jt 

where Jt is the thermal inertia of the thrust bearings, rt is the 
friction coefficient for the thrust bearings, and the Ht^i terms 
are heat transfer coefficients. The radial bearings behave sim- 
ilarly: 

Tr = — Hrp{Tr — To) — Hr,2{Tr ~ To)) 

where Jr is the thermal inertia of the radial bearings, is the 

friction coefficient for the radial bearings, and the terms 
are heat transfer coefficients. Note that rt and contribute 
to the overall friction coefficient r. 

The overall input vector u is given by 

U{t) = [ps{t) Pd{t) Ta{t) V{t) Ws(f)]^. 

The measurement vector y is given by 

y{t)=[uj{t) Q{t) Tt{t) Tr{t) To(f)]^. 

Fig. 5 shows nominal pump operation. The input voltage (and 
frequency) are varied to control the pump speed. The electro- 
magnetic torque is produced initially as slip is 1. This causes 
a rotation of the motor to match the rotation of the magnetic 
field, with a small amount of slip remaining, depending on 
how large the load torque is. As the pump rotates, fluid flow 
is created. The bearings heat up as the pump rotates and cool 
when the pump rotation slows. 

Damage Modeling 

The most significant forms of damage for pumps are impeller 
wear, caused by cavitation and erosion by the flow, and bear- 
ing failure, caused by friction-induced wear of the bearings. 
In each case, we map the damage to a particular parameter in 
the nominal model, and this parameter becomes a state vari- 
able in x(f) that evolves by a damage progression function. 
These functions are parameterized by a set of unknown wear 
parameters, forming the unknown parameter vector 9{t). 
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Time (hours) Time (hours) 

Radial Bearing Temperature Bearing Oil Temperature 



0 1 2 3 0 1 2 3 

Time (hours) Time (hours) 

Figure 5. Nominal pump operation. 

Impeller wear is represented as a decrease in impeller area 
A [7,9]. We use the erosive wear equation [12]. The ero- 
sive wear rate is proportional to fluid velocity times friction 
force. Fluid velocity is proportional to volumetric flow rate, 
and friction force is proportional to fluid velocity. We lump 
the proportionality constants into the wear coefficient wa to 
obtain 

A = -waQ^- 

A decrease in the impeller area will decrease the pump pres- 
sure, which, in turn, reduces the delivered flow, and, there- 
fore, pump efficiency. The pump must operate at a certain 
minimal efficiency. This requirement defines an EOT crite- 
ria. We define A~ as the minimum value of the impeller 
area at which this requirement is met, hence, Teol = 1 if 
A{t) < A~. 

Bearing wear is captured as an increase in friction. Sliding 
and rolling friction generate wear of material which increases 
the coefficient of friction [1,12]: 

rt(t) = WtTtUJ^ 
fr{t) = lUrTruA , 

where wt and Wr are the wear coefficients. The slip com- 
pensation provided by the electromagnetic torque generation 
masks small changes in friction, so it is only with very large 
increases that a change in lo will be observed. These changes 
can be observed much more readily through the bearing tem- 
peratures. Limits on the maximum values of these tempera- 
tures define EOL for bearing wear. We define and as 
the maximum permissible values of the friction coefficients. 


Algorithm 1 SIR Eilter 

Inputs: yfe 

Outputs: 

for i = 1 to TV do 

ei^p{ek\ei_A 

w’t ^p(yfc|x|,6>*j,,ufe) 

end for 

N 
i = l 

for i = 1 to TV do 
wi ^ w{/W 

end for 

^ Resample({(x^,ey,«;i,}[Li) 

before the temperature limits are exceeded over a typical us- 
age cycle. So, Teol = 1 if > r:^ or rr{t) > r+. 
Vibration and acceleration sensors have also been used in 
pumps for bearing monitoring, e.g., in [8], however, when 
using such methods it is difficult to map changes in vibra- 
tion back to changes in the thrust bearings, radial bearings, or 
both, while also quantifying the amount of damage. 

4. Damage Estimation 

In model-based prognostics, damage estimation reduces 
to joint state-parameter estimation, i.e., computation of 
p(xfc, 0fc|yo:fc). A general solution to this problem is the 
particle filter, which may be directly applied to nonlinear 
systems with non-Gaussian noise terms [13]. In particle fil- 
ters, the state distribution is approximated by a set of discrete 
weighted samples, called particles. 

With particle filters, the particle approximation to the state 
distribution is given by 

where N denotes the number of particles, and for particle i, 
denotes the state vector estimate, 6], denotes the parameter 
vector estimate, and w\. denotes the weight. The posterior 
density is approximated by 

TV 

p(xfc, 0 fc|yo:fc) « 0 ,^)(dxfcd 0 fc), 

i^l 

where 8(^^i^ Qi^'^{(b^kd6k) denotes the Dirac delta function lo- 
catedat 

We use the sampling importance resampling (SIR) particle fil- 
ter, using systematic resampling [14]. The pseudocode for a 
single step of the SIR filter is shown as Algorithm 1 . Each 
particle is propagated forward to time k by first sampling 
new parameter values, and then sampling new states using the 
model. The particle weight is assigned using y^. The weights 
are then normalized, followed by the resampling step [13]. 

Here, the parameters 9k evolve by some unknown process 
that is independent of the state x^. However, we need to 
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Figure 6. ^ adaptation scheme. 


assign some type of evolution to the parameters. The typi- 
cal solution is to use a random walk, i.e., for parameter 9, 
6k = 9k-i + where is sampled from some dis- 
tribution (e.g., zero-mean Gaussian). With this type of evolu- 
tion, the particles generated with parameter values closest to 
the true values should be assigned higher weight, thus allow- 
ing the particle filter to converge to the true values. 

The selected variance of the random walk noise determines 
both the rate of this convergence and the estimation perfor- 
mance once convergence is achieved. Therefore, it is very 
desirable to tune this parameter to obtain the best possible 
performance. A large random walk variance will yield quick 
convergence but tracking with too wide a variance, whereas 
too small a random walk variance will yield a very slow con- 
vergence, if at all, but, once achieved, tracking will proceed 
with a very small variance. One approach is to use kernel 
shrinkage, in which the random walk noise is diminished over 
time [15]. This approach assumes that the parameter is con- 
stant, but in reality, this may not be the case, so some amount 
of noise should still be included to account for unmodeled de- 
viations in the parameter value over time. In [16], this noise 
(viewed as a hyper-parameter) is tuned using outer correction 
loops based on prediction error. In this case, the underlying 
prognostic model is assumed to contain only a single fault 
dimension, therefore it cannot be applied in our case. 

We develop a ^ adaptation method similar to [16], but with 
some key distinguishing features. First, we consider a multi- 
dimensional damage space, therefore, we must simultane- 
ously adapt the random walk noise for multiple parameter 
values. Second, we cannot use prediction error to drive the 
adaptation, because we cannot, in general, map errors in pre- 
diction to specific wear parameters, since each output is de- 
pendent on multiple damage mechanisms. Instead, we try to 
control the variance of the hidden wear parameter estimate to 
a user-specified range by modifying the random walk noise 
variance. Since the random walk noise is artificial, we should 
reduce it as much as possible, because this uncertainty prop- 
agates into the EOL predictions. So, controlling this uncer- 
tainty helps to control the uncertainty of the EOL prediction. 
Reducing the variance of the wear parameter can reduce the 
variance of the EOL prediction by several factors, and the im- 
provement is substantial over long time horizons. 

The algorithm for the adaptation of the ^ vector is given as 


Algorithm 2 ^ Adaptation 


Inputs: 

State: a 
Outputs: 
if fc = 0 then 
a ^ 0 
end if 

for all j e {1, 2, . . . , n£i} do 
Vj ^ RMAD({0* (i)}[I,) 
if a( j ) = 0 and Vj < T then 

a(i) ^ 1 

end if 

if a( j ) = 0 then 


else 
end if 

end for 


GO 




Algorithm 2, and Eig. 6 shows how it interacts with the parti- 
cle filter. We assume that the ^ values are tuned initially based 
on the maximum expected wear rates, e.g., if the pump is ex- 
pected to fail no earlier than 100 hours, then this corresponds 
to particular maximum wear rate values. The initial wear rate 
estimate values may start at 0. We use the relative median 
absolute deviation (RMAD) as the measure of variance: 


RMAD(A) = 100 


Median^ {\Xi — Median^ ( Aj ) | ) 
Median^ (Xj ) 


where A is a data set and Xi is an element of that set. We 
use RMAD because it is statistically robust, and, since it is a 
relative measure of spread, it can be treated equally for any 
wear parameter value. The adaptation scheme resembles a 
proportional control law, where the error between the actual 
RMAD of a parameter 0{j), denoted as Vj in the algorithm, 
and the desired RMAD value (e.g., 10%), denoted as v* in the 
algorithm, is normalized by Vj . The error is then multiplied 
by a factor P (e.g., 1 x 10“^), and the corresponding variance 
4(j) is increased or decreased by that percentage. We utilize 
two different setpoints. Eirst, we allow for a convergence pe- 
riod, with setpoint v*q (e.g., 50%). Once Vj reaches T (e.g., 
1.2u*g), we mark it using the a(j) flag, and begin to control 
it to a new setpoint v*^ (e.g., 10%). 

Because there is some inertia to the process of Vj changing 
in response to a new value of ^(j), the gain P cannot be too 
large, otherwise Vj will not converge to the desired value, in- 
stead, it will continually shrink and expand. In our experi- 
ments, P = 1 X 10“^ worked well over the entire range of 
values considered for each wear parameter. Ideally, the wear 
parameter variance would be zero, but the particle filter needs 
some amount of noise to accurately track the parameter. So, 
V* cannot be too small, and we have found that controlling to 
an RMAD of 10% introduces an acceptable amount of uncer- 
tainty while allowing for accurate tracking. 
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Algorithm 3 EOL Prediction 


Inputs: 


Outputs: {EOL\, 
for i = 1 to Af do 

k <— kp 


lAT 
, Ji=l 


e 


ei 


k ^kp 

while TeOl{^1; 
Predict u^. 


0 


k + l 


^p{ek+i\ei) 


xLl ~p(xfe+i|x^,6»'j^,Ui;) 

k^k + 1 

j *v^ 

-Xfc ^ 
ai j u* 

^ ^k+1 

end while 

EOLi ^ k 

kp 

end for 


5. Prediction 

Prediction is initiated at a given time kp. Using the cur- 
rent joint state-parameter estimate, p{xkp , Okp |yo:fc_p ), which 
represents the most up-to-date knowledge of the system 
at time kp, the goal is to compute p{EOLkp\yo-.kp) and 
p{RULkp |yo:fej=)- As discussed in Section 4, the particle fil- 
ter computes 

N 

pixkp ,9kp\yo:kp) '^kp ,ei^ ) (dxkp dOkp ) . 

i=l 

We can approximate a prediction distribution n steps forward 
as [17] 

P{xkp+n: 9kpPn\yO:kp') ~ 

N 

'y ’^kpd(x.^^^^P'^^^^){dXkp+nddkp+n)- 
i=l 

So, for a particle i propagated n steps forward without new 
data, we may take its weight as Similarly, we can ap- 

proximate the EOL as 

N 

p{EOLkp\yo-.kp) -^wUpopi^^{dEOLkp). 

i=l 

To compute EOL, then, we propagate each particle forward 
to its own EOL and use that particle’s weight at kp for the 
weight of its EOL prediction. 

If an analytic solution exists for the prediction, this may be di- 
rectly used to obtain the prediction from the state-parameter 
distribution. An analytical solution is rarely available, so the 
general approach to solving the prediction problem is through 
simulation. Each particle is simulated forward to EOL to ob- 
tain the complete EOL distribution. The pseudocode for the 
prediction procedure is given as Algorithm 3 [1]. Each parti- 
cle i is propagated forward until Teol{x\,9\) evaluates to 
1; at this point EOL has been reached for this particle. 

Note that prediction requires hypothesizing future inputs of 
the system, u^, because damage progression is dependent on 



Figure 7. Simultaneous prediction of impeller wear and 
thrust bearing wear in the pump. The damage trajectories are 
coming out of the page, increasing in r*, decreasing in A, and 
increasing in t. 


the operational conditions. Eor example, in the pump, an in- 
creased rotation speed will cause bearing friction to increase 
at a faster rate, and will cause an increased pump flow, which, 
in turn, will cause impeller wear to increase at a faster rate. 
The choice of expected future inputs depends on the knowl- 
edge about operational settings and the type of information 
the user is interested in, e.g., for a worst-case scenario, one 
would consider the pump running at its maximum rotation. 

Eig. 7 shows results from the simultaneous prediction of im- 
peller wear and thrust bearing wear for N — 100 (not all 
trajectories are shown in the lower plot). Initially, the parti- 
cles have a very tight distribution of friction and impeller area 
damage values, but the distribution of the wear parameters, 
WA and Wrt, is relatively large. As a result, the individual 
trajectories are easily distinguishable as EOL is approached. 
Because the damage threshold is multi-dimensional, we show 
also the projections of the trajectories onto the damage-time 
planes. The projection onto the A-t plane (right) shows the 
progression of A towards the A~ threshold as a function of 
time. The projections stop when EOL is reached, and the ver- 
tical dotted lines connecting the projections to the time axis 
indicate individual EOL predictions. Similarly, the projec- 
tion onto the rt-t plane (bottom) shows the progression of rt 
towards the threshold as a function of time. The dotted 
lines connecting to the time axis indicate EOL predictions. 
Eor some particles, A~ is reached first, while for others, rf' 
is reached first. The different EOL values along with parti- 
cle weights form an EOL distribution approximated by the 
probability mass function shown in the upper plot. 
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6. Results 

In this section, we present simulation-based experiments to 
analyze the performance of the prognostics algorithm in the 
case of multiple damage progression paths. We first define the 
metrics used to evaluate the algorithm performance. We then 
provide detailed results for a single experiment to demon- 
strate the approach, followed by results summarized over a 
large number of experiments. 

Evaluation Metrics 

We evaluate the performance of the wear parameter estima- 
tion by quantifying estimation accuracy and spread. Accu- 
racy is calculated using the percentage root mean square error 
(PRMSE), which expresses relative estimation accuracy of w 
as a percentage: 


PRMSE^ = 100 

where Wk denotes the estimated wear parameter value at time 
k, denotes the true wear parameter value at k, and Mean^ 
denotes the mean over all values of k. In computing PRMSE, 
we ignore the initial time frame associated with convergence 
of the wear parameter estimate (from 0 hours up to 30% of 
the true EOL). 


Meanfc 

[ wk-wiy 


[v < J \ 


X 10 ■ 



j Mean(tU() 

Alin(tut) andMa^jut) 


0-P ^ ^ ^ ^ ^ ^ ^ 1 

0 5 10 15 20 25 30 35 40 


t (hours) 



Figure 8. Simultaneous estimation of pump wear parameters 
for N = 500, T = 60%, = 50%, = 10%, and P = 

1 X 10-3. 


We calculate the spread using RMAD as defined in Section 4. 
Eor estimation spread, for time k, we compute for wear pa- 
rameter w, RMADu, fc using the distribution of wear param- 
eter values given by the particle set at k as the data set. We 
denote the average RMAD over multiple k using: 

RMADu, = Meanfc(RMADu, fc). 

In computing estimation spread, we also ignore the initial 
time frame associated with convergence of the wear parame- 
ter estimate. 

For a particular prediction point kp, we compute measures of 
accuracy and spread for the prediction. For accuracy, we use 
the relative accuracy (RA) metric [10]: 

\RULl^-Metin,{RULl^)\ \ 

RULi^ ) ■ 

RA is averaged over each prediction point to obtain a single 
value that characterizes the overall accuracy, denoted as RA. 

We calculate prediction spread using RMAD, which we de- 
note as RMADflj/i for the RUE prediction. To obtain a sin- 
gle value for overall spread, RMAD is averaged over all pre- 
diction points starting from the prediction at which a prog- 
nostics horizon (where RA is within a specified bound) is 
first reached, denoted using RMAD/jf/L- Prognostics perfor- 
mance is summarized using the a-X metric which requires 
that for a given prediction time A, at least j3 of the RUE prob- 
ability mass lies within a of the true value [10]. 


Demonstration of Approach 

We first provide an example scenario to illustrate the ap- 
proach. Fig. 8 shows the estimation results for the hidden 
wear parameters, with w\ = 2 x IQ-^, wj" = 4 x and 

w* = 2x IQ-33. Initially, the estimate bounds are very large, 
however, as the estimates begin to converge, the RMAD of 
each is reduced to 50% through the adaptation scheme, and 
then to 10%. Once convergence has occurred, tracking pro- 
ceeds very well. The RMAD is maintained around 10% to the 
end of the experiment. The PRMSE of the different wear pa- 
rameters are correspondingly low, with PRMSE^,^ = 4.36 
PRMSE^, = 3.60, and PRMSE„^ = 5.51. The mean 
RMADs of the wear parameters are RMAD^,^ = 8.60, 
RMAD^j = 8.42, and RMAD„,^ = 8.29, which are less 
than the controlled value of 10%. 

Prediction performance is shown by the a-X plot of Fig. 9. 
Impeller wear damage dominates the EOF prediction. The 
accurate and precise wear parameter estimates yield corre- 
spondingly accurate and precise RUE predictions. Here, 
a = 0.1 and [3 = 0.5, so the a-X test requires that 50% 
of the probability mass lies within 10% of the true value at 
each prediction point. The test succeeds at all but the last pre- 
diction point, although the probability mass contained within 
the a-bounds, 49.6%, is very close to the requirement of 50%. 
The average RA is 97.16%. The average RMAD of the RUE 
distribution is 9.14%. Maintaining the variance of the wear 
parameter estimates maintains also the RMAD of the RUE 
(though not necessarily to the same setpoint). 
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Table 1. Estimation and Prediction Performance 


n 

PRMSEu,^ 

PRMSEs, 

PRMSE^^ 

RMADu,^ 

RMADujj 

RMADu,^ 

RA 

RMADfl;[/T 

1 

6.44 

6.64 

4.45 

8.44 

8.38 

8.30 

96.17 

10.24 

10 

5.38 

2.64 

3.25 

8.55 

8.76 

8.53 

96.79 

10.68 

100 

4.60 

2.71 

2.40 

9.12 

8.82 

8.88 

95.99 

11.65 
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Figure 9. a-A performance with a = 0.1 and f3 = 0.5 for 
N = 500, T = 60%, = 50%, = 10%, and P = 

1 X 10-3. 


Simulation Results 

We performed a number of simulation experiments in which 
combinations of wear parameter values were selected ran- 
domly within a range, with N = 500. We selected values in 
[0.5 X 10-3, 4 X 10-3] at increments of 0.5 x 10-3 for wa, in 
[0.5 X 10-^^, 7 X 10-^^] at increments of 0.5 x 10-^^ for wt, 
and in [0.5 x 10-^^, 7 x 10-^^] at increments of 0.5 x 10-^^ 
for Wr, such that the maximum wear rates corresponded to a 
minimum EOL of 20 hours. In order to conhrm that the wear 
parameter variance could still be maintained with additional 
sensor noise, we varied the sensor noise variance by factors of 
1, 10, and 100, and performed 20 experiments for each case. 
We considered the case where the future input of the pump is 
known, and it is always operated at a constant RPM. Hence, 
the only uncertainty present is that involved in the noise terms 
and that introduced by the particle hltering algorithm. 

The averaged estimation and prediction performance results 
are shown in Table 1. In all experiments, we used T = 60%, 
= 50%, = 10%, and P = 1 X 10-3. In each of 

the cases, the PRMSE for the different wear parameter esti- 
mates remained at most around 6.6% for the normal amount 
of noise, and under 5% for increased noise. We attribute the 
higher PRMSE of the normal noise cases to a couple out- 
lier scenarios where convergence was slower, throwing off 
the estimate early on. In these cases, the median PRMSEs 
were under 5%. The PRMSE for wa is on average higher 


than that for the bearing wear parameters because the flow 
measurement Q is relatively more noisy than the temperature 
measurements Tt and Tr- 

The RMAD of each wear parameter was successfully con- 
trolled to 10%, averaging around 8 to 9%. This trans- 
lated to good prediction performance, with the RA averaging 
around 96% and the RMAD of the RUE prediction averaging 
around 11%. Even as the noise increases, the variance con- 
trol scheme was able to maintain the RMAD setpoint, and so 
RMADi{( 7 L increased only slightly as sensor noise increased. 

Eig. 10 shows the RMAD of the wear parameters as a function 
of wear parameter value. Here, it is clear that the RMAD can 
be controlled well independently of the wear parameter value. 
Performance is similar across different wear parameters and 
their values, translating to the similar prediction performance 
observed across different wear parameter values. 

7. Conclusions 

We investigated the issues of multiple damage progression 
paths and developed a model-based prognostics methodology 
to accommodate them. Damage progression paths are char- 
acterized by a fault or damage variable and a set of wear pa- 
rameters that describe how they evolve in time. Particle biters 
perform joint state-parameter estimation in order to estimate 
the health state of the component. The state-parameter dis- 
tribution is then extrapolated to the EOL threshold to com- 
pute EOL and RUE predictions in the presence of multiple 
damage progression paths. A novel variance control mecha- 
nism keeps the uncertainty necessary for proper functioning 
of the particle biter in check, in order to maintain the uncer- 
tainty of the unknown wear parameters at a desired level. The 
framework was applied to a centtifugal pump, and the results 
demonstrated good performance over a range of wear param- 
eter values and sensor noise levels. 

In higher dimensional systems, the particle biter requires a 
very large number of particles to track successfully. Using 
only 500 particles was sufficient for good results here, but 
as the number of states or damage mechanisms needed to 
be tracked increases, the number of particles must increase 
also. Eor large N, the particle biter approach may not be 
efficient enough. In future work, we would like to investi- 
gate alternative approaches with reduced computational bur- 
den for high-dimensional state spaces. Also, the model-based 
approach presented here could possibly be complemented by 
data-driven methods that utilize pump vibration or accelera- 
tion sensors. 
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Figure 10. RMAD of the wear parameter as a function of 
wear parameter value. 
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